Skip to content

Add glm_pred fused GLM predictor virtual track - #87

Open
aviezerl wants to merge 38 commits into
masterfrom
feat/glm-pred
Open

Add glm_pred fused GLM predictor virtual track#87
aviezerl wants to merge 38 commits into
masterfrom
feat/glm-pred

Conversation

@aviezerl

Copy link
Copy Markdown
Collaborator

Summary

  • New virtual track type glm_pred that computes fused generalized linear model predictions at each genome position entirely in C++
  • R API: glm_pred.create(), glm_pred.info(), glm_pred.ls(), glm_pred.rm()
  • C++ engine with track grouping, scaling groups, selector-track per-bin model selection, per-bin weight skipping, and fast logistic transforms

Performance

Full mm10 genome (128.7M positions at 20bp iterator): 8.7 min with 63 cores on /dev/shm. 285x faster than equivalent R-level gextract with computed expressions.

Design

Each entry in the predictor follows: smooth → scale (cap + normalize) → transform (logistic) → weight. Interactions compute entry_i × entry_j / scale_factor, optionally transformed. A selector track bins each position to choose from K weight columns (for stratified models like per-GC-bin LASSO).

Key optimizations:

  • Track groups: read one track's super-window via mmap once, process all shifted sub-windows from the buffer
  • Scaling groups: aggregate+scale once per unique (track, shift), apply multiple transforms
  • Per-bin weight skipping: after reading the selector, skip track groups where all entries have zero weight for that bin
  • Fast exp(): polynomial approximation for logistic transforms (exact std::exp kept for lse aggregation)

Test plan

  • 171 tests covering: parameter validation, pipeline correctness (sum/lse), scaling, logistic transforms, multi-entry weighting, interactions, edge cases, sparse tracks, multi-chromosome, unaligned windows, selector tracks (K=2 basic, out-of-range NaN, interactions, per-bin bias, backward compat)
  • Full test suite: 18135 pass
  • Genome-wide inference validated against independent R-level reconstruction (r=1.000) and Tamar's pipeline (r=0.999)

@aviezerl
aviezerl force-pushed the feat/glm-pred branch 12 times, most recently from 1984137 to d795595 Compare April 21, 2026 10:42
aviezerl and others added 16 commits June 1, 2026 21:42
Adds a new "glm.predict" virtual track type that evaluates a fused
generalized linear model (per-GC-bin LASSO) at each genome position in
a single pass over the underlying motif and GC tracks.

- R API: glm_pred.create / .ls / .rm / .info plus a vignette and
  ~930 lines of tests covering scaling, logistic transforms, kernels,
  GC interactions, multi-bin prediction and sourceless evaluation.
- C++: GlmVarProcessor evaluates N motif entries (cap+normalize
  scaling, up to multi-parameter logistic transforms, kernel weights)
  plus optional tile-tile interactions, with bin-major weight layout
  for stride-1 cache access and shared track handles to avoid FD
  exhaustion under parallel eval.
- Wiring: GLM_PREDICT val_func in TrackExpressionVars,
  TrackExpressionScanner hooks, pkgdown reference section, and a
  FixedBin helper needed for kernel windows.
Single C++ call that replaces the R training-side pipeline
(gextract chunks -> scale_motif_energy -> pivot_wider -> logistic
transforms -> GC interactions) with one pass over the motif and GC
tracks, returning the complete n_peaks x n_features matrix directly.

Implementation (src/GlmFeatureExtractor.{h,cpp}):
- Tiles peaks internally, no 1.8M-row intermediate data frame.
- Opens all motif tracks + GC track once per chromosome, reads via
  mmap with float-precision lse_accumulate (bit-identical match with
  misha's vtrack LSE; two-pass double LSE gave ~4e-5 error amplified
  through the logistic chain).
- Supports both dense (FixedBin) and sparse tracks; sparse lookup
  uses binary search over the chromosome index instead of a linear
  scan (GC track has ~100M intervals per chromosome).
- Applies cap+normalize scaling and the four logistic heads in a
  single pass with no intermediate allocations.

R wrapper glm_extract_features() assigns column names for the motif,
GC and GC-interaction blocks. Comparison tests check exact match
against the R pipeline on real motif tracks.

For a 33K-peak / 191-motif workload this replaces ~10-20 min of
R-side processing with ~32s single-threaded.
Computes genome-wide quantiles (e.g. p=0.9999 for per-motif cap
values used by glm_extract_features) for a batch of motif tracks in
parallel without going through gquantiles / virtual tracks.

Implementation (src/GlmBatchQuantiles.cpp):
- One std::thread per track (default min(n_tracks, hw_concurrency,
  40)); each thread opens all chromosomes via mmap, computes LSE at
  each iterator position and collects values.
- Exact quantile via std::nth_element instead of
  StreamPercentiler's approximation.
- Avoids R/vtrack overhead and doMC fork setup; ~1.7x faster per
  track than gquantiles.

R wrapper glm_batch_quantiles() returns an n_tracks x n_percentiles
matrix. Comparison tests check values against gquantiles (diff
~0.02-0.06 at p=0.9999, within the approximate-vs-exact gap) and
are skipped when no local misha DB is available.
Phase 1 of batched multi-track functions plan. The old
GlmBatchQuantiles.cpp custom worker pool is replaced with the shared
BatchTrackScan<Reducer> skeleton; TopKQuantile reducer is defined in
BatchQuantiles.cpp and operates in fallback mode (full-vector storage
+ nth_element), matching pre-refactor behavior bit-for-bit.

- src/BatchTrackScan.{h,cpp,tpp}: shared templated scan driver with
  per-(track, chrom) task queue, sliding max/min deque (array-backed),
  monotonic interval-mask cursor, needs_pruning / needs_lower_bound
  constexprs, and strict no-R-API worker contract.
- src/BatchQuantiles.cpp: TopKQuantile reducer + C_gquantiles_multi
  .Call entry (renamed from C_glm_batch_quantiles).
- src/misha-init.cpp, R/misha-package.R, NAMESPACE: symbol rename.
- R/glm-features.R: dispatch updated to C_gquantiles_multi.

Existing test-glm-batch-quantiles.R passes (skipped in CI due to
missing mm10 db; R wrapper glm_batch_quantiles continues to work).
5 tests exercising the new BatchTrackScan code path on gdb.init_examples():
dense and sparse track paths, single and multiple percentiles, determinism
across thread counts. Complements test-glm-batch-quantiles.R which requires
the mm10 motif db and skips in CI.
- run_batch_scan now merges per-(track,chrom) task state into a per-track
  accumulator inside each worker (under a per-track mutex) and frees the
  task buffer. Prevents unbounded memory growth from holding all task
  buffers resident until main-thread merge. Returns BatchTrackScanResult
  (merged per_track_states + per-track error_messages) instead of raw
  tasks vector. Critical for Phase 1 fallback mode (full-value storage)
  and for Phase 2 median fallback on large workloads.
- Fix misleading "rewind" comment in BatchTrackScan.tpp — the code
  forward-jumps next_bin_to_push past a mask gap, not rewinds.
- Document sparse-path aggregator choices (double accumulator for
  SUM/AVG, no pruning).
- Document CAP invariant on SlidingExtremum::push_bin.

Tests: 15/15 pass on gdb.init_examples; 4 existing tests still skipped
on missing mm10 db.
TopKQuantile reducer now supports three modes:
  - use_fallback=true  -> full-vector storage + nth_element (Phase 1 behavior)
  - top-K heap mode    -> min-heap of top-K values (for p >= 0.5)
  - bottom-K heap mode -> max-heap of bottom-K values (for p < 0.5)
Heap built lazily when buf.size() reaches K. Merge across (track, chrom)
states concatenates and re-trims via nth_element — cheaper than N·K heap
pushes near K_MAX.

Adaptive K = ceil((1-min_p) N_est × 1.2), clamped to K_MAX=10M. Mixed-tail
percentiles and K>K_MAX trigger fallback with a warning.

Aggregator templating: func arg accepts lse|avg|sum|max|min, routed through
BatchTrackScan's WindowAggFunc enum and scan-driver dispatch. Sliding-max
(and optional sliding-min) pruning computes window upper/lower bounds and
calls Reducer::prune() before the aggregator math — skipping the exp/log
chain entirely for extreme quantiles.

Intervals support: optional _intervals data.frame; converted main-thread
into per-chrom sorted GInterval vectors; scan driver advances a monotonic
cursor to skip positions outside the mask. reducer.boundary() fires at gaps
(no-op for TopKQuantile; will matter for ThresholdScreen in Phase 4).

Signature: C_gquantiles_multi takes 9 args (added _func, _intervals). R
wrapper adds func="lse" and intervals=NULL defaults so existing callers
still work.

Tests: 12 new in test-batch-quantiles-phase2.R. Full suite 18,182 / 20 / 0.
Phase 2 review + new parity tests caught a real bug: top-K mode didn't
count pruned-but-valid positions toward n_total, so the rank calculation
in topk_finalize drifted from the fallback path. Concretely, at p=0.95
on gdb.init_examples dense_track, top-K produced 0.20667 while the
mixed-tail fallback produced 0.19333 — same scan, different rank.

Fix: State::count_pruned() increments n_total without pushing to the
heap. Driver calls it when prune() returns true AND window_max is
finite (at least one non-NaN bin → the aggregate would be non-NaN).

Also addresses review items:
- Rewrote misleading needs_lower_bound comment on TopKQuantile (actual
  value is true, original comment claimed false).
- Changed heap-build trigger from size==K to size>=K so a post-merge
  accumulator heapifies on its first accept rather than growing past K.
- Renamed "top-K clamp triggers fallback warning for p=0.5" test (the
  example db doesn't trigger the clamp; test actually asserts no warn).

Two new parity tests now assert bit-identical agreement between top-K
and fallback paths at the same p (via the mixed-tail fallback trick).
These would have caught the counting bug earlier.

Tests: 29/29 pass in batch-quantiles. Full suite 18,184 / 20 / 0.
Adds BatchSummary.cpp with a Summary reducer on BatchTrackScan<Summary>
and wires gsummary to dispatch:
  - Single expression (any R expression, including character scalars)
    always takes the legacy path — preserves bit-exact back-compat with
    the ~30 existing gsummary regression tests.
  - Character vector of length > 1 routes to .detect_fast_path, which
    accepts bare track names and simple vtracks (func in {lse, avg, sum,
    max, min}, single source, 1D iterator, consistent sshift/eshift).
    On dispatch it calls C_gsummary_multi and returns a data.frame with
    columns {track, n, n_nan, min, max, sum, mean, sd}.
  - Slow path for multi-expression arbitrary expressions is deferred to
    Phase 5; current behavior is to error with a clear message.

R/batch-dispatch.R hosts the dispatch helpers (.detect_fast_path,
.describe_single_expr, .fast_dispatch_msg) that will also be used by
gscreen (Phase 4) and gquantiles (Phase 6).

For bare tracks without an explicit iterator, the fast path uses the
track's native bin size as both iterator and window so semantics match
gsummary's implicit per-bin iteration.

Tests: test-batch-summary.R adds 8 cases (single-expr back-compat,
multi-track shape, parity with per-track slow-path, intervals mask,
vtrack-windowed, complex-expr fall-through, multi-expr slow-path error).
Full suite: 18,209 pass / 20 skip / 0 fail.
Adds BatchScreen.cpp with a ThresholdScreen reducer on
BatchTrackScan<ThresholdScreen>. Each task accumulates passing scan
positions into runs (merging consecutive positions within iterator-step
distance) and flushes at interval-mask gaps via boundary(). The run
flushes as [cur_start, cur_end + iterator_step), matching legacy
gscreen bin semantics (a passing position at pos emits a bin-width
interval [pos, pos + bin_size)).

Pruning: GT/GE skip when window_max < threshold; LT/LE skip when
window_min > threshold; EQ has no useful bound and always runs the
aggregate.

R/batch-dispatch.R gains .detect_screen_fast_path and .screen_op_to_int
helpers. The parser requires each expression to be a simple
"<lhs> <op> <const>" comparison where <op> in {<, <=, ==, >=, >} and
<lhs> satisfies .describe_single_expr (bare track or vtrack). All
expressions must share the same (func, sshift, eshift) tuple.

gscreen R wrapper extended: vector-of-expressions with length > 1 and
no intervals.set.out route to C_gscreen_multi. Single-expression calls
always take the legacy path to preserve back-compat. Multi-expression
slow path errors cleanly (Phase 5 scope). The returned data.frame is
long-form with columns {chrom, start, end, track}; the track column is
rewritten on the R side to show the original input expression strings
(not underlying source track names).

Tests: test-batch-screen.R adds 7 cases covering single-expr legacy
preservation, multi-track shape, per-track parity with legacy gscreen,
all 5 comparison operators, interval-mask boundary flush (no spurious
fusion across mask gaps), multi-expr slow-path error, and intervals.set.out
incompatibility with multi-expr.

Full suite: 18,237 pass / 20 skip / 0 fail.
gquantiles now accepts a character vector of expressions (length > 1)
when fast=TRUE, returning a data.frame with columns {track, <percentiles>}.
Single-expression calls with fast=TRUE route through the same
C_gquantiles_multi path (returns a named numeric vector, matching legacy
shape). Single-expression calls with fast=FALSE (the default) preserve
legacy behavior exactly.

Default is fast=FALSE for single expressions to avoid silently changing
numeric outputs for existing callers (StreamPercentiler vs exact
nth_element can differ by 0.02–0.06 at p≥0.999 per the design doc).
When a single-expression call WOULD be fast-path eligible, an opt-in
informational message fires once per session pointing the user at
fast=TRUE. The flip to fast=TRUE as default will come in a later
release, announced in NEWS.

Multi-expression without fast=TRUE errors cleanly with a message
pointing at the scope (Phase 5 slow-path multi-expr is out of scope
for this round).

Roxygen for gquantiles expanded with two detail sections documenting
(a) the numerical drift between fast and slow paths and (b) the memory/
speed behavior near p=0.5 (fallback to full storage).

Tests: test-gquantiles-dispatch.R adds 7 cases (single-expr default
back-compat, single-expr fast=TRUE shape, multi-expr fast=TRUE
data.frame, multi-expr without fast=TRUE error, multi-expr with complex
expr error, fast-vs-slow parity at matched iterator, vtrack vector).

Full suite: 18,257 pass / 20 skip / 0 fail.
Phase 4 review caught a real bug: when two or more input expressions
shared the same underlying source track (e.g. c("t > 1", "t > 2")),
the R-side track-column rewrite collapsed all rows onto the first
expression. The fallback was position-counting by consecutive runs of
matching track names, which never advanced because every row had the
same name string.

Fix: C_gscreen_multi now returns a track_idx integer column (0-based
into the input vector) in place of the track name. The R wrapper
indexes expr[track_idx + 1] — unambiguous regardless of duplicates
and simpler than the previous two-branch mapping.

Also addressed review I1 (dispatch overhead): added a cheap lexical
gate .looks_like_bare_name() to .detect_fast_path so legacy callers
with arithmetic expressions don't pay the gtrack.info / .gvtrack.get
resolution cost on every call.

Tests: test-batch-screen.R adds a regression test (two thresholds on
the same underlying dense_track) that fails before the fix and passes
after. Full suite: 18,260 pass / 20 skip / 0 fail.
Review caught a real bug: BatchSummary::n was never being incremented
for NaN-aggregate windows, and the n_nan branch inside accept() was
unreachable because the driver filtered NaN before calling accept. Net
result: df$n_nan was always 0, and df$n undercounted by the number of
NaN-aggregate windows — divergent from legacy gsummary semantics.

Fix: add a State::nan_seen() reducer hook. The driver now calls
state.nan_seen() when aggregate_window returns NaN (or when the
sparse-path window has zero non-NaN bins), and state.accept(val, pos)
otherwise — each reducer decides what to do:

  - Summary::nan_seen: ++n; ++n_nan (matches legacy num_bins /
    num_nan_bins).
  - TopKQuantile::nan_seen: no-op (NaN values have no rank).
  - ThresholdScreen::nan_seen: flush_cur() (NaN <op> threshold is false
    for every operator in legacy gscreen, which breaks any passing run).

Summary::accept is simplified: it's now only reached for finite values,
so the isnan branch is removed.

Tests: test-batch-summary.R parity test now also asserts n == n_total
and n_nan == n_nan_legacy with tolerance=0. Fails without the fix,
passes with it. Full suite: 18,262 pass / 20 skip / 0 fail.
Two small optimizations to the batched multi-track scan hot path,
selected from a nine-item review after benchmarking each against a
warm-cache 5-track mm10 workload. Only the two that moved measurable
metrics (or were algorithmically obvious) are retained here.

1. Sparse-path monotonic cursor. `scan_sparse_inner` was doing a
   per-position binary search into `intervals` to find the first
   overlapping sparse interval. Since the outer `c` advances
   monotonically, `win_s` does too, so the search front can only move
   forward — replaced with a sticky `sparse_cursor` advance. Not
   measurable in the dense-track benchmark, but algorithmically clear
   and matches the pattern already used for the interval-mask cursor.

2. Runtime gate on sliding-deque maintenance. TopKQuantile fallback
   mode (mixed-tail / K_MAX-clamped) never prunes, yet the driver was
   still pushing every bin through the sliding max/min deques every
   position. Added `State::pruning_active()` hook; hoisted the check
   once per task into a local bool; wrapped the entire deque-and-prune
   block in `if (pruning_enabled)`. ThresholdScreen returns false for
   EQ predicates (EQ never prunes).

   Measured B_fallback (mixed-tail quantile): 8.85s → 6.66–7.08s across
   stable runs (-20 to -25%). Small-interval screen/quantile cases
   regress ~30–40ms (+10%), but those cases are dominated by per-track
   NFS file-open latency so the absolute impact is minor.

   The runtime check is hoisted via a constexpr-branched IIFE so that
   for reducers with `needs_pruning=false` (Summary) the `if constexpr`
   still elides the block entirely.

Tests: full suite passes (one unrelated multicontig test-environment
flake under parallel runner, passes standalone).
Symmetric to the existing `prune()` hook. When the aggregator bound
already guarantees the threshold predicate passes (e.g. `window_min >
threshold` for GT with LSE/AVG/SUM/MIN), the driver skips
aggregate_window entirely and calls `accept_certain_pass(pos)` — which
extends/starts the passing run without re-deriving a value (threshold
screens only emit intervals, not values).

New reducer interface:
  - static constexpr bool supports_certain_pass
  - State::certain_pass(float upper, float lower)
  - State::accept_certain_pass(int64_t pos)

Only ThresholdScreen opts in; TopKQuantile and Summary need the actual
aggregated value so they set supports_certain_pass=false and the driver
elides the branch via if constexpr.

Pruning and certain-pass are mutually exclusive at any given position
(prune requires bound fails predicate; certain-pass requires bound
satisfies predicate), so the order of checks inside the driver loop
doesn't matter for correctness.

Note: the fast path only fires for aggregators whose lower-bound (or
upper-bound for LT/LE) is informative. MAX's lower bound is -inf and
MIN's upper bound is +inf (see aggregate_lower_bound / aggregate_upper_
bound), so certain_pass on `vtrack_max > t` or `vtrack_min < t` never
fires — that's correct: those aggregators can only fail via prune, not
pass via certain_pass.

Benchmark (5 vtrack LSE screens on mm10, threshold -1e6 so every
window passes):
  - before: 37.6s
  - after:   3.65s  (~10× speedup)

Correctness: output of fast-path vs legacy slow-path is bit-identical
(266 passing intervals, same chrom/start/end). Full test suite passes
(18,262 tests).
aviezerl added 21 commits June 1, 2026 21:42
Pre-existing issues introduced by the phase-3/4/6 commits that only
surfaced now that CI re-ran R-CMD-check on a PR refresh.

1. .tpp file triggers "unlikely file names for src files" WARNING.
   R CMD check only recognizes a fixed set of source extensions; .tpp
   isn't one of them. Renamed src/BatchTrackScan.tpp →
   src/BatchTrackScan_impl.h (still included by BatchTrackScan.h);
   updated include-guard and one comment cross-reference.

2. Codoc mismatches for four Rd files (Rd was missing args that exist
   in the code):
    - gquantiles: Rd missing `fast`. (@param already present in source
      roxygen but Rd was stale.)
    - gsummary:   added `@param fast` to roxygen.
    - gscreen:    added `@param fast` to roxygen.
    - glm_batch_quantiles: added `@param func` and `@param intervals`.

Ran alutil::style_and_document() to regenerate all Rd files and apply
styler. Full batch test suite passes locally.
…sitional order

Chromids were computed in R via match() against gintervals.chrom_sizes(intervals),
which only reflects the ordering of chromosomes present in the input subset, not
the chromkey insertion order. When the input intervals lacked any chromosome that
comes earlier in the chromkey (e.g. chrM/chrY in mm10), every chromosome after
the missing one was silently shifted by one — chrX became id 19 = chrM, returning
all-zero features because chrM is far shorter than chrX peak coordinates.

Pass chromosome names through to C++ and resolve them with chromkey.chrom2id()
inside C_glm_extract_features. Per-chrom calls happened to be unaffected (a
single-chrom subset's positional id always matched the chromkey id), so the
regression test asserts that multi-chrom output equals concatenated per-chrom
output and that chrX rows are non-zero.
Add an n_threads argument to glm_extract_features (default
getOption("gmax.processes", 1L); 0 = auto: min(n_peaks, hardware_concurrency, 40)).

Per-chrom strategy: open all motif + GC track handles once on the main
thread, pre-materialize sparse track intervals/values, then fan out the
peak loop across worker threads via an atomic chunk counter (CHUNK=256).
Output rows are disjoint per peak so there is no write contention, and
shared track handles are read-only after materialization.

Worker threads must not longjmp out via verror() — open_track_static and
the worker lambda throw TGLException/std::runtime_error instead, which the
main thread catches and re-raises with verror() once all threads have
joined. aggregate_lse / aggregate_sum become static so they can be called
without holding the extractor instance.

Adds bit-identity test across n_threads ∈ {1, 4, auto} on 600 multi-chrom
peaks. C_glm_extract_features arity bumps 13 → 14 in misha-init.cpp.
Prepare for multi-selector support. Replace scalar glm_selector_* fields
with std::vector<...> forms of length M. Field names and semantics will be
exercised by the parsing and runtime updates in subsequent commits.
Read 'selector_tracks' as a character vector (length M >= 1) and
'selector_breaks' as a list of M numeric break vectors. Build M
BinFinders and compute column-major strides (first selector varies
fastest). Validate Pi K_m matches glm_num_bins set by the weights matrix.
start_chrom now iterates over var.glm_selector_track_names (length M)
and populates glm_selector_handles/fixedbins/bin_sizes vectors.
For M=0 (no stratification) the loop is a no-op.
Read M selector mmap pointers per position, run BinFinder per selector,
and assemble the compound bin via sum_m b_m * stride_m. Any selector
producing NaN, out-of-range, or a missing pointer propagates NaN to the
output (any-failure -> NaN semantics, matching the prior single-selector
contract).

Local variable named 'nsel' (not 'M') to avoid colliding with the
existing 'int M = (int)var.glm_interactions.size();' declared a few
lines later in the same scope.
Replace selector_track (singular) and selector_breaks (numeric vector)
with selector_tracks (character vector, length M) and selector_breaks
(list of M numeric vectors). K_total is the product of per-selector
K_m. Existing single-selector callers must migrate to a 1-element
selector_tracks and a 1-element selector_breaks list.
Two roxygen lines and one stop() error message still mentioned the
singular 'selector_track' after the API rename. Updated to
'selector_tracks' for consistency.
Mechanical rewrite of all call sites and assertions to use
selector_tracks (character vector) and selector_breaks (list of
numeric vectors). Validation-message expectations updated:
'at least 2' -> 'length >= 2' to match the new wording in
R/glm-pred.R. The 'selector_breaks = 0.5' scalar test is wrapped
as 'selector_breaks = list(0.5)' to exercise the per-element
length check rather than the new is.list(...) check that fires
first.
Three new test_that blocks plus a small .val2bin_finder helper that
emulates BinFinder::val2bin(right=TRUE, include_lowest=TRUE) via
cut(..., include.lowest=TRUE, right=TRUE). The helper avoids the
findInterval-vs-BinFinder boundary mismatch I hit when first writing
the planned tests.

(1) Cartesian product across two selectors with K_total = 6 strata,
    hand-computed expected = bias[k+1] + alpha[1, k+1] * x where
    k = (b1 - 1) + (b2 - 1) * K1 (column-major flatten).
(2) Any-NaN / any-OOR selector value -> NaN output, with a positive
    assertion that some positions actually trigger the OOR path.
(3) M=1 multi-selector path matches a hand-computed single-selector
    reference.
Replace the single 'selector_track / selector_breaks' section with a
multi-selector formulation. Adds the column-major compound-bin formula
b(p) = sum_m (b_m(p) - 1) * prod_{m' < m} K_{m'}, an R snippet showing
array-to-matrix flatten via R's default column-major, and a one-line
note that the M=1 case is a length-1 selector_tracks plus length-1
selector_breaks list.
Multi-selector glm_pred.create() now takes weights/bias/interaction_weights
as labeled multi-axis arrays with dim = c(N, K_1, ..., K_M) (and
c(M_int, K_1, ..., K_M) for interactions). Trailing-axis names and bin
labels are auto-derived from selector_tracks and selector_breaks via
cut(include.lowest=TRUE, right=TRUE), so glm_pred.info() round-trips a
fully labeled object that you can index by stratum tuple at debug time.

This kills the "did I flatten in the wrong order?" bug class on the
pre-release multi-selector path. The flat N x prod(K_m) matrix shape is
no longer accepted for M >= 2; the validation error names the expected
shape and the migration. Single-selector (M = 1) still accepts a plain
matrix(N, K_1) and numeric(K_1) since there is no flatten ambiguity.
Scalar bias (default 0) is recycled to all strata.

C++ unchanged: bin-major [b*N + offset] indexing matches R's column-major
array storage of dim = c(N, K_1, ..., K_M).
Drop the sshift < eshift requirement so a shifts entry can also
translate the interval (sshift == eshift) instead of only expand it.
length(selector_tracks) >= 1 and length(br) >= 2 are validated above,
so K_per[m] >= 1 for every m and prod(K_per) >= 1 always.
Commit 3c9d354 dropped the sshift < eshift requirement so a shifts
entry can also translate the interval (sshift == eshift) and reverse
orientation. The validation test was not updated and has been failing
on branch HEAD ever since. Replace with positive tests that pin the new
semantics: pure-translation shift c(50, 50) and reverse-orientation
shift c(100, -100) should both create vtracks without error.
Five latent bugs found by a deep review pass on feat/glm-pred. Each fix
has at least one regression test (TDD-verified by reverting the fix and
watching the test fail).

* glm_pred.create now validates length of simple_cap,
  interaction_trans_family, and inter_trans_params (length 1 or N,
  matching the entry-side trans_family policy). Wrong lengths previously
  passed silently: simple_cap entries beyond the input length got the
  default FALSE; interaction_trans_family beyond input got NA via R
  out-of-bounds subset, silently disabling the transform on those
  interactions.

* glm_pred virtual tracks now verror on entry-track open exceptions
  (TrackExpressionVars::start_chrom). Previously the catch silently
  nulled cached_fixedbin/cached_sparse, masking real corruption or
  removed-track errors as all-NaN output for the affected chromosome.
  The legitimate "no data on this chrom" path returns nullptr from
  create_and_init_1d_track without throwing and is still handled
  silently in the else-branch above the catch.

* glm_extract_features derives transform column-name suffixes from
  names(transforms). The default transforms is now a named list, so
  default usage produces the same column names as before. The previous
  magic L-value fingerprint silently relabeled non-canonical transform
  sets as t1..tN.

* gsummary multi-track fast path reads C output columns by name instead
  of position. C side already names the columns; this removes the silent
  column-swap risk if column ordering ever changes C-side. New
  test-batch-summary.R contract test pins the C/R column-name agreement.
With breaks like c(-Inf, x, Inf) every adjacent diff is Inf, the
equality check in BinFinder::init passes, and m_binsize stays Inf.
val2bin's uniform-binsize fast path then evaluates Inf/Inf = NaN;
(int)NaN is undefined behaviour and on x86_64 yields INT_MIN, which
the surrounding min(INT_MIN - 1, K - 1) wraps to K - 1. Result: every
value silently routes to the last bin -- glm_pred selectors with such
breaks (and any other BinFinder consumer) misclassify silently.

Specifically observed: a glm_pred vtrack with a c(-Inf, 0, Inf)
polarization selector ignored the negative-polarization slice of the
weights array; only the positive slice was ever read.

Force the binary-search path when m_binsize is non-finite. Added
regression test exercising c(-Inf, x, Inf) on a glm_pred selector;
test fails without the fix and passes with it. Broader vtrack-glm-pred
+ gpartition/gquantiles/gscreen/gsynth suites stay green.
The batch fast path used by gquantiles(fast=TRUE), gsummary(fast=TRUE),
and gscreen(fast=TRUE) computed the per-position window as
[c + sshift, c + eshift] anchored at a single iterator-step start c.
The slow path -- per Iterator_modifier1D::transform in
TrackExpressionVars.h -- spans the iterator interval [s, e] and applies
shifts to its boundaries: [s + sshift, e + eshift] = [c + sshift,
c + iterator_step + eshift]. The fast-path window was therefore
iterator_step shorter than intended.

Most-visible failure: a vtrack with gvtrack.iterator(sshift=0, eshift=0)
collapses the fast-path window to [c, c]; the `if (win_e <= win_s)
continue;` guard then skips every position, gquantiles returns NaN,
gsummary reports n=0/all-NaN, gscreen returns 0 hits. With non-zero
shifts the fast path returned silently-wrong values (window 20bp short
of the slow path's). Bare tracks with iterator != bin.size were also
wrong; the helper happened to mask this when iterator == bin.size by
returning eshift = bin.size, which only worked for that aligned case.

Fix:

- src/BatchTrackScan_impl.h: change win_e = c + eshift to
  win_e = c + iterator_step + eshift in both scan_fixedbin_inner and
  scan_sparse_inner. Widen n_bins_per_window to a safe upper bound
  (ceil(W/B) + 1 where W = iterator_step + eshift - sshift) so the
  LSE/SUM upper-bound pruning const stays correct for windows unaligned
  to the bin grid -- otherwise ThresholdScreen prunes valid passes (the
  symptom was gscreen returning 12 vs slow-path's 7825 even after the
  win_e fix on its own).

- R/batch-dispatch.R: .describe_single_expr for bare tracks now returns
  sshift = eshift = 0 with a separate default_iterator = bin_size
  (used only when iterator = NULL). Vtracks set default_iterator = NA
  to force explicit iterator. .detect_fast_path and
  .detect_screen_fast_path consume the new field. Updates the contract
  doc to spell out window semantics.

Tests: added vtrack-default-shifts, vtrack-non-zero-shifts, and
bare-track-iterator-not-bin-size cases to test-gquantiles-dispatch.R,
plus a vtrack-default-shifts gscreen case to test-batch-screen.R. The
existing "matches legacy on small examples" test still passes (it used
iterator == bin_size, the aligned case where the old code happened to
be correct). Full test suite green (18532 pass / 0 fail).

No NEWS entry: the fast=TRUE feature itself was introduced on this same
unreleased feat/glm-pred branch, so this bug never reached a release.
Defer the version bump until release; branch stays at 5.9.1 (master's
current release) with glm_pred / glm_extract_features / glm_batch_quantiles
and the multi-expression gsummary/gscreen/gquantiles fast paths recorded
under a development-version NEWS section.
A review of the batched multi-track fast path (gscreen/gsummary/gquantiles
fast=TRUE) for fast-path-vs-slow-path divergence turned up six issues. All are
on this unreleased branch, so no NEWS entry. Each is covered by a new parity
test (test-batch-fast-parity.R) that compares the fast path against the legacy
path on controlled tracks (signed values, NaN gaps, non-aligned chrom sizes)
and was confirmed red before / green after.

C++ (src/BatchTrackScan.cpp, BatchTrackScan_impl.h, BatchScreen.cpp):

1. SUM aggregate bounds were not valid bounds. aggregate_upper/lower_bound
   returned window_max/window_min * n_bins, but n_bins (= ceil(W/B)+1) is an
   over-estimate of the true non-NaN bin count and the per-bin extreme can have
   the wrong sign. This mis-pruned/certain-passed gscreen("...sum </> t") on
   signed tracks (missed AND spurious intervals) and skewed gquantiles(func=sum)
   pruning, even on all-positive tracks via the lower bound. Now
   max(window_max,0)*n / min(window_min,0)*n, which are universally valid.

2. All-NaN windows could certain-pass / prune-fuse instead of breaking the run.
   A window with no non-finite... no non-NaN bins has wmax == -inf; for MIN>t
   (upper bound +inf) and MAX<t (lower bound -inf) the non-informative bound
   landed on the predicate side and certain_pass fired, fusing the two data
   runs around a NaN gap into one interval. The driver now skips prune /
   certain-pass when wmax == -inf and falls through to aggregate_window ->
   nan_seen, matching the slow path (flush on NaN).

5. gscreen flushed run ends as cur_end + iterator_step with no clamp, so on a
   chromosome whose size is not a multiple of the iterator the last interval
   ran past the chromosome end. Clamp to chrom size at emit (the slow path
   clips the final iterator interval to the chrom end).

7. The batch path reads raw mmap bins and only tested isnan, while the legacy
   read path converts +/-inf to NaN. aggregate_window, the sliding deques, and
   the sparse path now treat non-finite (isfinite==false) bins as missing, so
   tracks containing inf bins agree with the slow path.

8. The sliding deque advanced its front AFTER pushing the new window's bins,
   momentarily holding ~2x the window across a step transition, which could
   alias past WINDOW_CAP for large iterators. Advance the front before pushing;
   the final deque state is identical but the live span is bounded to the
   window size.

R (R/batch-dispatch.R):

3. A bare sparse track was fast-path eligible with iterator=NULL, scanning a
   1bp grid instead of the irregular sparse iterator the slow path uses.
   Sparse bare tracks now require an explicit iterator (default_iterator = NA).

4. An explicit intervals scope was accepted regardless of grid alignment, but
   the fast path scans a global (origin-0) grid with point membership and only
   matches the slow path's clipped iterator when each scope interval is aligned
   to the iterator step. .scope_grid_aligned() now declines non-aligned scopes
   (allowing an end == chrom size, which both paths clip identically).

Not addressed here: gquantiles fast=TRUE uses nearest-rank (nth_element) rather
than the slow path's linear interpolation between order statistics. That is
documented behavior, not a bug; flagged separately for a decision on whether to
add interpolation for parity.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant